
# Sanghyun Hong, February 23, 2019

#rm(list = ls())
###################################
## Simulation Environment Setup  ##
## Simulation Environment Setup  ##
#################################################################
#pathname<-"C:/Temp/SimulationProject/"
effectSize_List<-c(0.00, 0.03, 0.06);
sigH_List<-c(0.006);
PubBias_List<-c(0.5);
MetaStudyN_List<- c(40,60,80,100);
simN<-simN;
#################################################################
#
#
#
#
#
#################################################################
# Import Packages and set seed
#################################################################
library(metafor) # rma; ramdon effect estimator
library(sqldf) # This pakcage is for SQL
set.seed(3);
#################################################################
#
#
#
#
#
#################################################################
# Load Packages and Set Path
#################################################################
library(splines)
library(stats)
library(ggplot2)
library(gridExtra)
library(cowplot)
pathnameAK<-paste(pathname,"SimCodes/AK_Codes/",sep="")
setwd(pathnameAK)
source("Modified_AllApplications.R")
#################################################################
#
#
#
#
#
###############################################
## Primary Study Data (logOR) Generation     ##
## Primary Study Data (logOR) Generation     ##
#################################################################
PrimaryStudy_LogOR <- function(TreatmentEffect, sigH, bias, primaryObs){
  if(bias==0){
    P <- 0.1 + TreatmentEffect + rnorm(1, 0, sigH);
  }else{
    P <- 0.1 + TreatmentEffect;
  }
  nT<-primaryObs;
  mT<-sum(rbinom(nT, size=1, prob=P));
  nC<-primaryObs;
  mC<-sum(rbinom(nC, size=1, prob=0.1)); 
  
  correction<-(mT*mC==0)*10^(-5);
  LogOR<-log(((mT+correction)/(nT-mT+correction))/((mC+correction)/(nC-mC+correction)));
  LogORse<-sqrt( (1/(mT+correction)) + (1/(nT-mT+correction)) + (1/(mC+correction)) + (1/(nC-mC+correction)) )
  LogORci<-LogOR+LogORse*qt(c(0.025, 0.975), df=10^100)
  LogORt<-LogOR/LogORse
  return(c(TreatmentEffect,sigH, round(log(((.1+TreatmentEffect)/(.9-TreatmentEffect) )/(.1/.9)),2), LogOR, LogORse, LogORci, (LogORci[1]>0), primaryObs))
}
###################################
## Meta Analysis Data Generation ##
## Meta Analysis Data Generation ##
#################################################################
MetaStudy <- function(effect, sigH, m, bias, PrimaryN){
  MetaStudyData<-matrix(0,nrow=m, ncol=10)
  colnames(MetaStudyData)<-c('StudyID','effect','sigH','TrueLogOR','EstimatedLogOR','StdErrLogOR','LowerBound','UpperBound','Pos_SigEffect', 'PrimaryStudyOBS')
  biasedStudy<-sample(c(ceiling(m*bias), floor(m*bias)),1)
  for(i in 1:m){
    primaryObs<-PrimaryN[i];
    if(runif(1, 0, 1)<bias){
      PrimaryStudy<-PrimaryStudy_LogOR(effect, sigH, bias, primaryObs)
      while(PrimaryStudy[8]==0){PrimaryStudy<-PrimaryStudy_LogOR(effect, sigH, bias, primaryObs)}
      MetaStudyData[i,]<-c(i,PrimaryStudy)
    }else{
      PrimaryStudy<-PrimaryStudy_LogOR(effect, sigH, bias, primaryObs)
      MetaStudyData[i,]<-c(i,PrimaryStudy)
    }}
  return(MetaStudyData)
}
#################################################################
#################################################################
#
#
#
#
#
###################################
## Creating Tables for Output    ##
## Creating Tables for Output    ##
#################################################################
matrixRow<-length(effectSize_List)*length(sigH_List)*length(MetaStudyN_List)*length(PubBias_List)*simN
FE_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(FE_Output)<-c('TrueLogOR','m','sigH','PubBias','EstimatedLogOR','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
RE_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(RE_Output)<-c('TrueLogOR','m','sigH','PubBias','EstimatedLogOR','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
WAAP_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(WAAP_Output)<-c('TrueLogOR','m','sigH','PubBias','EstimatedLogOR','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
WAAP2_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(WAAP2_Output)<-c('TrueLogOR','m','sigH','PubBias','EstimatedLogOR','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
PETPEESE_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(PETPEESE_Output)<-c('TrueLogOR','m','sigH','PubBias','EstimatedLogOR','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
AK_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(AK_Output)<-c('TrueLogOR','m','sigH','PubBias','EstimatedLogOR','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
AK3_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=16))
colnames(AK3_Output)<-c('TrueLogOR','m','sigH','PubBias','EstimatedLogOR','stdErr','tstat','pvalue','lowerBound','upperBound','tau2','H2','I2','bias','coverage','reportingBias')
SubTable<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=17))
colnames(SubTable)<-c('TrueLogOR','m','sigH','FppCnt','WappCnt','WappTotal','WappSample','WappFrac','Wapp2Cnt','Wapp2Total','Wapp2Sample','Wapp2Frac','AK1P','AK3P1','AK3P2','AK3P3','I2')
AveBias_Output<-as.data.frame(matrix(NA, nrow=matrixRow, ncol=6))
colnames(AveBias_Output)<- c('TrueLogOR','m','sigH','PubBias','EstimatedLogOR','Bias')
#################################################################
#################################################################
#
#
#
#
#
###################################
## Simulation Begins Here        ##
## Simulation Begins Here        ##
################################################################
start.time <- Sys.time();
sim_i<-0;
PrimaryN<-rep(c(50,100,100,250,500),100)
for(PubBias in PubBias_List){
  for(m in MetaStudyN_List){
    for(effectSize in effectSize_List){
      for(sigH in sigH_List){
        cnt<-0;
        for(simnulationLoop in 1:simN){
          cnt<-cnt+1;
          sim_i<-sim_i+1;
          data<-as.data.frame(MetaStudy(effectSize, sigH, m, PubBias, PrimaryN));
          ##################
          ## Average Bias ##
          ## Average Bias ##
          #########################################################
          AveBias_Output[sim_i,]<-c(data$TrueLogOR[1],m,sigH,PubBias,mean(data$EstimatedLogOR), mean(data$EstimatedLogOR)-data$TrueLogOR[1])
          #########################################################
          #########################################################
          #
          #
          #
          #
          #
          #
          #
          ###################
          ## FIXED EFFECTS ##
          ## FIXED EFFECTS ##
          #########################################################
          effectFE<-as.matrix(data$EstimatedLogOR)
          wFE<-1/data$StdErrLogOR^2
          RegFE<-summary(lm(effectFE~1, weights=wFE))
          u_FE <- as.numeric(RegFE$coefficients[1])
          seFE <- as.numeric(RegFE$coefficients[2])
          p_valueFE <- as.numeric(RegFE$coefficients[4])
          ciFE <- as.numeric(confint(lm(effectFE~1, weights=wFE), level=0.95))
          covFE <- (ciFE[1]<=data$TrueLogOR[1])*(ciFE[2]>=data$TrueLogOR[1])
          FE_Output[sim_i,]<-c(data$TrueLogOR[1],m,sigH,PubBias,u_FE,seFE,u_FE/seFE,p_valueFE, ciFE[1], ciFE[2], 0, 0, 0, u_FE-data$TrueLogOR[1], covFE, mean(effectFE-data$TrueLogOR[1]))
          rm(effectFE, wFE, RegFE, u_FE, seFE, p_valueFE, ciFE, covFE)
          #########################################################
          #########################################################
          #
          #
          #
          #
          #
          #
          #
          ###################
          ## Random Effect ##
          ## Random Effect ##
          #########################################################
          effectRE<-as.matrix(data$EstimatedLogOR)
          wRE<-1/data$StdErrLogOR^2
          RegRE<-summary(rma.uni(data$EstimatedLogOR~1,se=data$StdErrLogOR, intercept=TRUE, data=data, weighted=TRUE, method="DL", level=95, digits=5));
          tau2RE <- as.numeric(RegRE$tau2)
          I2RE <- as.numeric(RegRE$I2)
          H2RE <- as.numeric(RegRE$H2)
          u_RE <- as.numeric(RegRE$b)
          seRE <- as.numeric(RegRE$se)
          p_valueRE <- as.numeric(RegRE$pval)
          ciRE <- as.numeric(c(RegRE$ci.lb, RegRE$ci.ub))
          covRE<-(ciRE[1]<=data$TrueLogOR[1])*(ciRE[2]>=data$TrueLogOR[1])
          RE_Output[sim_i,]<-c(data$TrueLogOR[1],m,sigH,PubBias,u_RE,seRE,u_RE/seRE,p_valueRE, ciRE[1], ciRE[2], tau2RE, H2RE, I2RE, u_RE-data$TrueLogOR[1], covRE, mean(effectRE)-data$TrueLogOR[1])
          rm(effectRE, wRE, RegRE, H2RE, tau2RE, u_RE, seRE, p_valueRE, ciRE, covRE)
          #########################################################
          #
          #
          #
          #
          #
          #
          #
          ###################
          ## WAAP          ##
          ## WAAP          ##
          #########################################################
          wappcnt<-1;
          u_WAAP <- as.numeric(summary(lm(data$EstimatedLogOR~1, weights=1/(data$StdErrLogOR^2)))$coefficients[1])
         
          APSdata<-subset(data, (abs(u_WAAP/2.8)>=data$StdErrLogOR))
          if(nrow(APSdata)<2){APSdata<-data; wappcnt<-0;}
          APSdata<-cbind(APSdata, u_WAAP=u_WAAP)

          wapptot <- nrow(data);
          wappsam <- nrow(APSdata);
          wappfrac <- wappsam/wapptot;
          
          effectAPS<-as.matrix(APSdata$EstimatedLogOR)
          wAPS<-1/APSdata$StdErrLogOR^2
          RegAPS<-summary(lm(effectAPS~1, weights=wAPS))
          u_APS <- as.numeric(RegAPS$coefficients[1])
          seAPS <- as.numeric(RegAPS$coefficients[2])
          p_valueAPS <- as.numeric(RegAPS$coefficients[4])
          ciAPS <- as.numeric(confint(lm(effectAPS~1, weights=wAPS), level=0.95))
          covAPS <- (ciAPS[1]<=data$TrueLogOR[1])*(ciAPS[2]>=data$TrueLogOR[1])
          if(nrow(data)==nrow(APSdata)){APS_Bias<-u_APS-data$TrueLogOR[1];}else{APS_Bias<-u_APS-data$TrueLogOR[1];}
          WAAP_Output[sim_i,]<-c(data$TrueLogOR[1],m,sigH,PubBias,u_APS,seAPS,u_APS/seAPS,p_valueAPS, ciAPS[1], ciAPS[2], 0, 0, 0, APS_Bias, covAPS, APS_Bias)
          
          rm(u_WAAP, APSdata, effectAPS, wAPS, RegAPS, u_APS, seAPS, p_valueAPS, ciAPS, covAPS, APS_Bias)
          #########################################################
          #
          #
          #
          #
          #
          #
          #
          ###################
          ## WAAP2         ##
          ## WAAP2         ##
          #########################################################
          wapp2cnt<-1;
          u_WAAP <- as.numeric(summary(lm(data$EstimatedLogOR~1+data$StdErrLogOR, weights=1/(data$StdErrLogOR^2)))$coefficients[1])
          
          APSdata<-subset(data, (abs(u_WAAP/2.8)>=data$StdErrLogOR))
          if(nrow(APSdata)<2){APSdata<-data; wapp2cnt<-0;}
          APSdata<-cbind(APSdata, u_WAAP=u_WAAP)
          
          wapp2tot <- nrow(data);
          wapp2sam <- nrow(APSdata);
          wapp2frac <- wapp2sam/wapp2tot;
          
          effectAPS<-as.matrix(APSdata$EstimatedLogOR)
          wAPS<-1/APSdata$StdErrLogOR^2
          RegAPS<-summary(lm(effectAPS~1, weights=wAPS))
          u_APS <- as.numeric(RegAPS$coefficients[1])
          seAPS <- as.numeric(RegAPS$coefficients[2])
          p_valueAPS <- as.numeric(RegAPS$coefficients[4])
          ciAPS <- as.numeric(confint(lm(effectAPS~1, weights=wAPS), level=0.95))
          covAPS <- (ciAPS[1]<=data$TrueLogOR[1])*(ciAPS[2]>=data$TrueLogOR[1])
          if(nrow(data)==nrow(APSdata)){APS_Bias<-u_APS-data$TrueLogOR[1];}else{APS_Bias<-u_APS-data$TrueLogOR[1];}
          WAAP2_Output[sim_i,]<-c(data$TrueLogOR[1],m,sigH,PubBias,u_APS,seAPS,u_APS/seAPS,p_valueAPS, ciAPS[1], ciAPS[2], 0, 0, 0, APS_Bias, covAPS, APS_Bias)
          rm(u_WAAP, APSdata, effectAPS, wAPS, RegAPS, u_APS, seAPS, p_valueAPS, ciAPS, covAPS, APS_Bias)
          #########################################################
          #
          #
          #
          #
          #
          #
          #
          ###################
          ## PET-PEESE     ##
          ## PET-PEESE     ##
          #########################################################
          fppcnt <- 0
          effectFPP<-as.matrix(data$EstimatedLogOR)
          wFPP<-1/data$StdErrLogOR^2
          se1FPP<-data$StdErrLogOR
          se2FPP<-data$StdErrLogOR^2
          RegFPP<-summary(lm(effectFPP~1+se1FPP, weights=wFPP))
          ciFPP <- as.numeric(confint(lm(effectFPP~1 + se1FPP, weights=wFPP), level=0.95)[1,1:2])
          if(as.numeric(RegFPP$coefficients[1,3])>1.96){
            fppcnt <- 1
            RegFPP<-summary(lm(effectFPP~1+se2FPP, weights=wFPP))
            ciFPP <- as.numeric(confint(lm(effectFPP~1 + se2FPP, weights=wFPP), level=0.95)[1,1:2])
          }
          
          u_FPP <- as.numeric(RegFPP$coefficients[1,1])
          seFPP <- as.numeric(RegFPP$coefficients[1,2])
          p_valueFPP <- as.numeric(RegFPP$coefficients[1,4])
          covFPP <- (ciFPP[1]<=data$TrueLogOR[1])*(ciFPP[2]>=data$TrueLogOR[1])
          PETPEESE_Output[sim_i,]<-c(data$TrueLogOR[1],m,sigH,PubBias,u_FPP[1],seFPP[1],u_FPP[1]/seFPP[1],p_valueFPP[1], ciFPP[1], ciFPP[2], 0, 0, 0, u_FPP[1]-data$TrueLogOR[1], covFPP,mean(effectFPP)-data$TrueLogOR[1])          
          
          rm(effectFPP, wFPP,se1FPP, se2FPP, RegFPP, ciFPP, u_FPP, seFPP, p_valueFPP, covFPP)
          #########################################################
          #
          #
          #
          #
          #
          #
          #
          ###################
          ## A-K           ##
          ## A-K           ##
          #########################################################
          AKdata<-as.data.frame(cbind(data$EstimatedLogOR, data$StdErrLogOR, data$StudyID))
          attach(setup(AKdata, TRUE), warn.conflicts = FALSE)
          source("EstimatingSelection.R")
          p_valueAK <- 2*pt(abs(Psihat[1]/se_robust[1]), df=(nrow(my.data)-1),  lower.tail= FALSE);
          ciAK<-as.numeric(Psihat[1]+qt(c(.025, .975), df=(nrow(my.data)-1))*se_robust[1])
          covAK<-(ciAK[1]<=data$TrueLogOR[1])*(ciAK[2]>=data$TrueLogOR[1])
          
          AK_Output[sim_i,]<-c(data$TrueLogOR[1],m,sigH,PubBias,Psihat[1],se_robust[1],Psihat[1]/se_robust[1], p_valueAK, ciAK[1], ciAK[2], 0, 0, 0, Psihat[1]-data$TrueLogOR[1], covAK, mean(data$EstimatedLogOR)-data$TrueLogOR[1])
          Ak1p1<-Psihat[3]
          rm(Psihat, se_robust, p_valueAK, covAK, ciAK)
          detach("setup(AKdata, TRUE)")
          #########################################################
          #
          #
          #
          #
          #
          #
          #
          ###################
          ## A-K 3         ##
          ## A-K 3         ##
          #########################################################
          AKdata<-as.data.frame(cbind(data$EstimatedLogOR, data$StdErrLogOR, data$StudyID))
          attach(setup(AKdata, FALSE), warn.conflicts = FALSE)
          source("EstimatingSelection.R")
          p_valueAK <- 2*pt(abs(Psihat[1]/se_robust[1]), df=(nrow(my.data)-1),  lower.tail= FALSE);
          ciAK<-as.numeric(Psihat[1]+qt(c(.025, .975), df=(nrow(my.data)-1))*se_robust[1])
          covAK<-(ciAK[1]<=data$TrueLogOR[1])*(ciAK[2]>=data$TrueLogOR[1])
          
          AK3_Output[sim_i,]<-c(data$TrueLogOR[1],m,sigH,PubBias,Psihat[1],se_robust[1],Psihat[1]/se_robust[1], p_valueAK, ciAK[1], ciAK[2], 0, 0, 0, Psihat[1]-data$TrueLogOR[1], covAK, mean(data$EstimatedLogOR)-data$TrueLogOR[1])
          Ak3p1<-Psihat[3]; Ak3p2<-Psihat[4]; Ak3p3<-Psihat[5];
          rm(Psihat, se_robust, p_valueAK, covAK, ciAK)
          detach("setup(AKdata, FALSE)")
          #########################################################
          SubTable[sim_i,]<-c(data$TrueLogOR[1],m,sigH, fppcnt, wappcnt, wapptot, wappsam, wappfrac, wapp2cnt, wapp2tot, wapp2sam, wapp2frac, Ak1p1, Ak3p1, Ak3p2, Ak3p3, I2RE)
          rm(data)
          cat("Simulation Environment #1: ", sim_i, "/", matrixRow,"\n");
        }
 }}}}
#################################################################
#################################################################
#
#
#
#
#
###############################################
## Making Graphs and Saving Output           ##
## Making Graphs and Saving Output           ##
#################################################################
FE_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueLogOR, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from FE_Output group by TrueLogOR, m, sigH'))
RE_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueLogOR, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from RE_Output group by TrueLogOR, m, sigH'))
WAAP_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueLogOR, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from WAAP_Output group by TrueLogOR, m, sigH'))
WAAP2_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueLogOR, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from WAAP2_Output group by TrueLogOR, m, sigH'))
PETPEESE_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueLogOR, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from PETPEESE_Output group by TrueLogOR, m, sigH'))
AK_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueLogOR, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from AK_Output group by TrueLogOR, m, sigH'))
AK3_OUTPUT_FINAL <- as.data.frame(sqldf('select TrueLogOR, m, sigH, PubBias, avg(bias) as Bias, avg(bias*bias) as MSE, avg(tau2) as tau2, avg(H2) as H2, avg(I2) as I2, avg(coverage) as coverage, avg(reportingBias) as reportingBias from AK3_Output group by TrueLogOR, m, sigH'))

BIAS_FINAL <- cbind(FE_OUTPUT_FINAL$TrueLogOR, FE_OUTPUT_FINAL$m, FE_OUTPUT_FINAL$sigH, FE_OUTPUT_FINAL$PubBias,
                    FE_OUTPUT_FINAL$Bias, RE_OUTPUT_FINAL$Bias,
                    WAAP_OUTPUT_FINAL$Bias, WAAP2_OUTPUT_FINAL$Bias, PETPEESE_OUTPUT_FINAL$Bias, AK_OUTPUT_FINAL$Bias, AK3_OUTPUT_FINAL$Bias)
colnames(BIAS_FINAL)<-c("TrueLogOR","m","sigH","PubBias","FE_Bias","RE_Bias","WAAP_Bias","WAAP2_Bias","PETPEESE_Bias","AK_Bias","AK3_Bias")


MSE_FINAL <- cbind(FE_OUTPUT_FINAL$TrueLogOR, FE_OUTPUT_FINAL$m, FE_OUTPUT_FINAL$sigH, FE_OUTPUT_FINAL$PubBias,
                   FE_OUTPUT_FINAL$MSE, RE_OUTPUT_FINAL$MSE,
                   WAAP_OUTPUT_FINAL$MSE, WAAP2_OUTPUT_FINAL$MSE, PETPEESE_OUTPUT_FINAL$MSE, AK_OUTPUT_FINAL$MSE, AK3_OUTPUT_FINAL$MSE)
colnames(MSE_FINAL)<-c("TrueLogOR","m","sigH","PubBias","FE_MSE","RE_MSE","WAAP_MSE","WAAP2_MSE","PETPEESE_MSE","AK_MSE","AK3_MSE")


COV_FINAL <- cbind(FE_OUTPUT_FINAL$TrueLogOR, FE_OUTPUT_FINAL$m, FE_OUTPUT_FINAL$sigH, FE_OUTPUT_FINAL$PubBias,
                   FE_OUTPUT_FINAL$coverage, RE_OUTPUT_FINAL$coverage,
                   WAAP_OUTPUT_FINAL$coverage, WAAP2_OUTPUT_FINAL$coverage, PETPEESE_OUTPUT_FINAL$coverage, AK_OUTPUT_FINAL$coverage, AK3_OUTPUT_FINAL$coverage)
colnames(COV_FINAL)<-c("TrueLogOR","m","sigH","PubBias","FE_COV","RE_COV","WAAP_COV","WAAP2_COV","PETPEESE_COV","AK_COV","AK3_COV")


SubTable_Ave<-sqldf('select TrueLogOR, m, sigH, avg(FppCnt) as FppCnt, avg(WappCnt) as WappCnt, avg(WappTotal) as WappTotal, avg(WappSample) as WappSample, avg(WappFrac) as WappFrac, avg(Wapp2Cnt) as Wapp2Cnt, avg(Wapp2Total) as Wapp2Total, avg(Wapp2Sample) as Wapp2Sample, avg(Wapp2Frac) as Wapp2Frac, avg(AK1P) as AK1P, avg(AK3P1) as AK3P1, avg(AK3P2) as AK3P2, avg(AK3P3) as AK3P3, avg(I2) as I2 from SubTable group by TrueLogOR, m, sigH')

write.csv(BIAS_FINAL, paste(pathname,"/Output/SimulationEnvironment01_BIAS.csv", sep=""))
write.csv(MSE_FINAL,  paste(pathname,"/Output/SimulationEnvironment01_MSE.csv", sep=""))
write.csv(COV_FINAL,  paste(pathname,"/Output/SimulationEnvironment01_COV.csv", sep=""))
write.csv(SubTable_Ave, paste(pathname,"/Output/SimulationEnvironment01_Miscellaneous.csv", sep=""))

end.time <- Sys.time();
time.taken <- end.time - start.time;
cat("Total Time Taken: ", time.taken, "\n\n\n")
#################################################################
#################################################################